During the last 25 years, the oil and gas industry has sought to gain more subsurface property information from seismic data—in particular, detailed information about subsurface pore fluids, porosity, lithology, pressure and geometry. This has been accomplished, in part, by using a combination of seismic data and well log information to increase the accuracy of subsurface properties estimated from seismic data. These efforts have been rewarded by improved exploration success rates and development well “sweet spot” selection. This success is, in large part, based upon the recovery of relative acoustic impedances (product of density and compressional velocity changes at an interface).
FIG. 1 shows a number of subterranean formations 10 for which detailed subsurface information is sought. Subterranean formations 10 may include one or more reservoirs or producing zones. A seismic survey is taken over this subterranean region of interest using seismic sound sources 12 and receivers 14. The receivers 14 record sound waves which travel in wavefronts from sources 12 to the subterranean formations 10 where a portion of the sound waves are reflected back to receivers 14 and recorded as seismic traces. Often a recording vehicle 16 is used to store recorded data. Also, preferably, a wellbore 20 passes through the subterranean formations 10 and logging tools are run in the wellbore 20 to obtain subsurface property information about the elastic properties of the subsurface formations 10.
FIG. 2 is an exemplary schematic of a single shot profile including a single source 12 and multiple receivers 14 laid out in a 2D geometry. The schematic shows the ray path geometry (normal to the propagating waterfront) that primary (direct path) seismic reflection energy travels in going from a source location to reflection points 22 and then back to receiver locations. An angle θr is shown, which is the reflection angle for the largest source-receiver distance.
Shot profiles are the basic seismic survey components required to form a common mid-point (CMP) gathers as shown in FIG. 3. In this figure, the reflection angle θr is shown for the nearest source-receiver pair. A CMP trace gather is a collection of shot profile traces that have differing shot to receiver distances and a common surface or mid-point location. This location corresponds to a common subsurface reflection point 24 if geologic layers are not dipping. Because geologic layers are often complexly structured, CMP gather traces are processed using various imaging and noise suppression technologies into common reflection point (CAP) trace gathers. These gathers are the typical collection of traces that are input to amplitude inversion algorithms. The gathers contain information about how primary reflection strength or amplitude changes at a common geologic interface as a function of the ray path arrival angle. Many arrival angle traces are required to invert for the underlying interface properties that have given rise to the observed reflection amplitude.
FIG. 4 depicts a convolutional model of an amplitude versus offset (AVO) primary reflection signal and shows how a primary reflection signal amplitude of a trace can be related to interfaces 26 between geologic formations 10. Data from well logs and cores taken from wellbores 20 demonstrate that the physical properties of geologic formations typically differ from one formation to another. Sound waves propagating through these formations are most sensitive to a formation's elastic properties including compressional velocity, vp, shear velocity vs, and density properties ρ and, to a lesser degree, their absorptive properties. Differences in elastic properties at the interface between two formations govern the amplitude of the reflected and transmitted wavefront's amplitude relative to that of the downgoing or arriving wavefront's amplitude. As FIG. 4 indicates, the magnitude of reflected signal amplitude can be approximated by a weighted sum of three elastic reflectivities or differences in interface elastic properties normalized by an average of the interface's elastic properties. The weighting function for each reflectivity is a trigonometric function of the reflection angle as described by equation (1).A(θ)≅R0+G*sin2 74 +RP*sin2 θ tan2 θ  (1)where: A(θ) amplitude reflected at angle θ                θ [averaged interface angle]=(θrefl+θrefr)/2        R0 [p-impedance reflectivity]=ΔvP/2vP+Δρ/2ρ        G [gradient reflectivity]=ΔvP/2vP−K(Δρ/2ρ+ΔvS/vS);        K=(2vS/vP)2; and        RP [p-reflectivity]=ΔvP/2vP.        
To a propagating seismic wavefront, a geologic section made up of many formation interfaces will act like a series of reflectivities spaced apart by the time it takes for the sound wave to travel between the interfaces. The primary reflection signal response is calculated by replacing each of the reflectivities by a copy of the propagating wavelet, or wavefront disturbance, scaled by the magnitude and sign of the reflectivity. The surface recorded primary response is the sum of all the time shifted, scaled wavelets, or convolution of the wavelets, with the interface reflectivity series. The two arrows at the bottom of FIG. 4 indicate that the goal of AVO amplitude inversion is to convert AVO traces into their component reflectivities while the goal of the AVO attribute analysis is to imply changes in geologic properties at the interface that have generated the inverted reflectivities.
Well log data recorded in wellbores 20 drilled into geologic formations 10 can be used to estimate the elastic and absorptive properties of the geologic formations 10 and subsequently the primary reflection signal amplitude that is part of the total recorded seismic survey response. Well log data can be used to statistically characterize the expected reflection response and to provide processing quality measures as explained below.
The groundwork for angle dependent amplitude inversion was laid in the 1950's when Bortfeld, R., 1961, Approximations to the reflection and transmission coefficients of plane longitudinal and transverse waves: Geophys. Prosp., v. 9, p. 485-502, described a linearized expression for the Zoeppritz reflection coefficient equation. Lindseth, R. O., 1979, Synthetic sonic logs—a process for stratigraphic interpretation: Geophysics, 44, p 3-26, implemented the inversion of trace amplitude for impedance. A qualitative AVO analysis was implemented shortly thereafter by Ostrander, W. J., 1984, Plane-wave reflection coefficients for gas sands at nonnormal angles of incidence: Geophysics 49, 1637-1648.
Subsequently, both qualitative and quantitative amplitude versus offset (AVO) inversions have been employed to estimate subsurface geologic properties. Amplitude inversion comprises the process of predicting one or more of the component reflectivities (convolved with a known wavelet) giving rise to an AVO reflection response from a collection of common subsurface reflection point traces ordered by increasing reflection angle or shot to receiver offset.
FIG. 5 shows a single interface between two media 30 and 32, a reflection angle θrefl and refraction angles θrefr for an arriving impulsive wavelet having unit amplitude arriving and the sum and difference elastic properties that are important in determining the strength of the reflected amplitude pulse. Also shown is that media has properties vp1, vs1 and ρ1, and media 32 has properties vp2, vs2 and ρ2.
Table 1 provides the factors for determining reflection amplitude at an interface.
TABLE 1Factors for Determining Reflection Amplitude at An Interfaceinterface reflectivities = Δf(vS, vp, ρ)/g(vS, vp, ρ)vP = media compressional wave velocityvS = media shear wave velocityρ = media densityΔ vS = vS2 − vS1; vS = (vS2 + vS1)/2ΔvP = vP2 − vP1; vP = (vP2 + vP1)/2Δρ = ρ2 − ρ1; ρ = (ρ2 + ρ1)/2
Table 2 provides definitions for reflectivities that are commonly used in AVO inversion. Note that Δf( . . . ) and g( . . . ) can be different functions of the media's differential and averaged properties.
TABLE 2Definitions of Reflectivities Used in AVO Amplitude InversionReflectivitySymbolDefinitionp-impedanceR0ΔvP/2vP + ρΔ/2ρ = Δ(ρvP)/ρvPp-velocityRPΔvP/2vPhybrid shearRSH2(vS/vP)2 (ρΔ/2ρ + Δ vS/vS) − Δ(ρvS2)/(ρvP2)gradientGΔ vP/2 vP − (2vS/vP)2(ρΔ/2ρ + Δ vS/vS) =RP − 2RSHdensityRρΔρ/2ρs-impedanceRSIΔ vS/2 vS + ρΔ/2ρ = Δ(ρvS)/ρvSs-velocityRSΔ vS/2 vS
FIG. 6 is a plot of a typical AVO response for a single interface using the full Zoeppritz plane wave reflection coefficient Eqn. (1). The figure shows that for typical geologic interfaces and at small reflection angles, the reflection coefficient is nearly constant. At angles approaching 30°, the reflection coefficient typically diminishes a few percent and at angles approaching the critical angle, the reflection coefficient approaches unity. There are many approximations to the exact Zoeppritz equation. These approximations stem from the complex structure of the Zoeppritz equation and the types of angle dependent amplitude information available from the seismic data acquisition and processing methods.
The Aki-Richards equation, which is a linearized version of the plane wave Zoeppritz equation, yields an excellent approximation to primary reflection amplitude for normal to precritical incidence angles and for small contrast interfaces. It is also an equation upon which many amplitude inversion algorithms are based. In terms of the geometry shown in FIG. 5 and elastic media parameters, the equation can be written as:A(θ)≅R0+G*sin2 θ+RP*sin2 θ tan2 θ  (1)where: A(θ)=amplitude reflected at angle θ                θ [averaged interface angle]=(θrefl+θrefr)/2        θ [p-impedance reflectivity]=ΔVP/2VP+Δρ/2ρ        G [gradient reflectivity]=ΔVP/2VP−K(Δρ/2ρ+ΔVS/VS); K=(2VS/VP)2         RP [p-reflectivity]=ΔVP/2VP and VP, VS, ρ, ΔVP, ΔVS and Δρ are defined in Table 1.        
Eqn. (1) provides that AVO amplitude, as a function of subsurface reflection angle, is a sum of products of trigonometric functions of reflection angle and three reflectivities R0, G and RP. These reflectivities, in turn, are functions of differences and averages of compressional wave velocity, shear wave velocity and density properties across the subsurface interface. The gradient term G is a particularly complicated reflectivity involving the sum and products of three other reflectivities. In addition, the angle in Eqn. (1) is a function of the average of the incidence and refraction angles at the interface. In general Eqn. (1) can be formulated using other trigonometric angle functions and reflectivities having the form shown below:reflectivity=Δf(VS, VP, ρ)/g(VS, VP, ρ)  (2)
The measured amplitudes can be inverted to yield the three interface reflectivities by making reflection amplitude measurements for three or more source to receiver offsets and accurately measuring the overburden velocity field from offset dependent travel times to estimate reflection angle.
TABLE 3                    ⁢                  REFLECTIVITY        ⁢                                  ⁢        MODEL            _                          ⁢                  A        i            =                        R          0                +                              R            1                    ⁢                      SIN            2                    ⁢                      ∅            i                          +                              R            2                    ⁢                      SIN            2                    ⁢                      ∅            i                    ⁢                      TAN            2                    ⁢                      ∅            i                                          WHERE      :                          ⁢                          ⁢              A        i              =          MEASURED      ⁢                          ⁢      AMPLITUDE      ⁢                          ⁢      AT      ⁢                          ⁢              ∅        i                                ⁢                  R        0            =              CHANGE        ⁢                                  ⁢        IN        ⁢                                  ⁢        IMPEDANCE        ⁢                                  ⁢                  (                                    ΔV              P                        ⁢            ρ                    )                                        ⁢                  R        1            =              GRADIENT        ⁢                                  ⁢        COEFFICIENT                                ⁢                  R        2            =              CHANGE        ⁢                                  ⁢        IN        ⁢                                  ⁢                  V          P                                        ⁢                  Ø        i            =              REFLECTION        ⁢                                  ⁢        ANGLE                                ⁢                                        REFLECTOR            ⁢                                                  ⁢            SOLUTION                    _                ⁢                                  [                                                            R                0                                                                                        R                1                                                                                        R                2                                                    ]            =                                    [                                                            ANGLE                                                                              GEOMETRY                                                                                                  (                                          3                      ×                      3                                        )                                                                        ]                                -            1                          ⁡                  [                                                                      ΣA                  i                                                                                                                          ΣA                    i                                    ⁢                                      SIN                    2                                    ⁢                                      ∅                    i                                                                                                                                            ΣA                    i                                    ⁢                                      SIN                    2                                    ⁢                                      ∅                    i                                    ⁢                                      TAN                    2                                    ⁢                                      ∅                    i                                                                                ]                    
Table 3 depicts an unconstrained, least squares method (L2 norm) for amplitude inversion that is typically used to solve for reflectivities for a reflection whose amplitudes and angles have been estimated by data processing methods. The L2 unconstrained method minimizes an error function of the form:E2=Σ(Ai−(R0+R1*sin2 θi+R2*sin2 θi tan2 θi))2   (3)by solving: ∂E2/∂R0=0                ∂E2/∂R1=0        ∂E2/∂R2=0where: Ai and θi=event amplitudes and reflection angles.        
A goal of quantitative amplitude inversion is for the signal to noise ratio of inverted attributes to be comparable over small spatial distances to that of the stack section response. As those practiced in the art of quantitative seismic amplitude inversion are well aware, it is very difficult to obtain quantitatively useful amplitude inversion results with Eqn. (1) alone. Major sources of error typically include amplitudes contaminated with various noises and inaccurate estimates of incidence angle at large reflection angles. In noisy data with a limited reflection angle range and uncorrected signal distortion effects, an unconstrained amplitude inversion will generate inverted reflectivities that have very low S/N ratios compared to the stack section or p-impedance section.
In order to improve the S/N ratio of amplitude inversion results, Eqn. (3) has been formulated using other error norms (eg. an L1 norm) and various constraints that supplement the seismic amplitude information. The constraints can be “hard” constraints—those that change the form or the variables in the reflectivity equation or “soft” constraints—additional information that is included as part of the error function which is to be minimized.
An AVO equation used in the 1970's for shorter offset acquisition geometries and non-amplitude preserving processing sequences, modeled amplitude with offset with an equation of the form:A(X)≅R0+G*X2  (4)                where: X=shot to receiver offset.        
This formulation uses a “hard” constraint, relative to equation (1) that implies that the sin2 θ tan2 term is negligible and that within a single scalar constant and at small reflection angles, squared offset is a good proxy for the sin2 θ term. Unmigrated CMP gathers with NMO removed were the usual input to a least squares fit of R0 and G to processed amplitude. A time-averaged energy or envelope difference of the R0 and the G terms was used as a direct hydrocarbon indicator and no other constraint data were employed.
Later the reflection angle was calculated from a measured velocity field using a straight ray approximation resulting in:A(θst)≅R0+G*sin2 θst  (5)where: θst=straight ray approximation to the incidence angle=tan−1(X/(t0Vrms)).
Eqn. (5) implies that the sin2 θ tan2 θ term of equation (1) is negligible and that the straight ray angle is a good approximation for the subsurface reflection angle—a good assumption when there is little velocity acceleration. The use of this equation also initiated an analysis method for the detection of anomalous event behavior by crossplotting of an event's intercept against its gradient.
A form of Eqn. (1) that uses minimal hard constraint assumptions but requires amplitude information over an angle range spanning normal to critical angles and that has been used in generating the inversion examples used in this patent is:A(θ)≅R0−2RSH*sin2 θ+RP*tan2 θ  (6)                where: RSH=Δ(ρVS2)/ρVP2        
In order to compensate for inadequacies in the AVO data due to noise and distortion, assumptions regarding the relationships between a rock's compressional velocity, shear velocity and density and the form of the VP/VS term have led to versions of Eqn. (1) that reflect these “hard” constraints by altering the variables in the reflection equation. These hard constraints include empirical relationships like A+BVP==VS implying that:ΔVS/VS=(BVP/VS)ΔVP/VP  (7)                where A and B are constants.        
The Gardner rule governing the relationship between velocity and densityρ=C VpK  (8)where C and K are constants has been used to imply thatΔρ/ρ=K ΔVP/VP  (9).
Other formulations of Eqn. (1) can be expressed in terms of p-impedance and shear impedance reflectivity result when the parameter K=(2VS/VP)2 is set to a constant.
Soft constraints in the form of data weighting and damping constraints are also used in an amplitude inversion to account for noise and distortion in the seismic amplitude data. Soft constraints added to equation (3) lead to an error function of the form:E2=ΣWi(Ai−(R0+R1*sin2 θi+R2*sin2 θi tan2 θi))2+WCOR02+WC1R12+WC2R22+Wff(R0, R1, R2)2+  (10)where: Wi=signal to noise ratios estimates of amplitudes Ai;                WC0,C1,C2,f=damping factors applied to reflectivities; and        f(R0, R1, R2)=damping factors for sums or differences of reflectivities.        
A shortcoming associated with utilizing equation (10) is establishing criteria for choosing the weights and damping factors. Ideally damping weights should be zero and amplitude weighting factors equal to unity if the amplitude data have no noise or distortion. Because field recorded shot profiles often are more noise than signal and because the shot profile signal suffers from various distortion effects, the type of data processing described below is highly desirable prior to an amplitude inversion
A typical field recorded shot profile of seismic data consists of primary reflection signals significantly distorted by the acquisition system, the effects of transmission thru the earth from shot to receiver and the influence of shot generated, ambient and acquisition noises. FIGS. 7A & 7B are actual shot profiles after normal moveout (NMO) correction (correction for the velocity of sound wave propagation) displayed with (FIG. 7A) and without (FIG. 7B) an initial noise suppression step in which primary reflections should appear as parallel flat lying coherent events. As shown by FIGS. 7A & 7B, real data primary signal reflections are difficult to detect because of various additive noises (in this case ground roll and multiples) and signal distortion effects which would yield very noisy inverted reflectivities.
FIG. 8 is a schematic depicting some of the important factors that can distort a propagating seismic wavelet and its subsurface reflection behavior. Influencing factors include the seismic source strength, the source and receiver directivity, coupling, array characteristics and near surface layering, wavefront spreading or divergence losses, inelastic absorption, interbed multiple scattering, local geologic dip, the curvature of reflecting interfaces lateral earth heterogeneity, transmission coefficient losses, the recording system filters and the subsurface reflection coefficients which are to be determined.
The purpose of seismic data processing sequence, as shown below in Table 4, is to remove multiples and other noise from the seismic data and to compensate for the effects of acquisition and earth filters on primary reflections by applying various correction algorithms to the seismic survey data. Some of the steps may themselves consist of a sub-sequence of steps while other steps may be repeated more than once in the sequence with different parameters as the signal and noise structure of the data becomes more evident.
TABLE 4Processing/Inversion Sequence ComponentsPost-migration and pre-inversionGeneric Processing Stepsamplitude conditioning stepsacquisition geometry assignmentprestack imaged crp gathersinitial velocity analysisCRP noise suppression *data regularizationwavelet spectral equalization *shot/cmp noise suppression *residual event alignment *signal corrections earth/recordingspatial amplitude balancewavelet equalization/deconvolution *offset/angle amplitude balancefinal velocity analysisangle/offset selection mute *CMP noise suppression *amplitude inversion/uncertainty *prestack migration/velocityrefinement* = step may use well log statistics
An objective of an amplitude processing sequence (one that precedes an amplitude inversion) is to make various corrections to the data while also preserving the relative amplitude behavior of primary reflections in space and reflection angle or offset. This particular objective can also be quantitatively stated as that of recovering amplitudes within a single scalar constant of the earth's plane wave, band-limited subsurface reflectivity so that the response at each interface corresponds to the convolution of a wavelet having known relative amplitude, phase and timing with the local interface reflection coefficients.
An amplitude processing sequence may have 20 or more steps including those listed in Table 4. The left column of Table 4 includes typical steps in a generic processing sequence that produces migrated CRP gathers while the right column lists steps in a post-migration, pre-inversion amplitude conditioning sequence. Important generic processing steps include geometry assignment, velocity analysis, data regularization, passes of noise suppression in different domains, signal distortion corrections, a wavelet equalization correction, a final velocity analysis and a migration or imaging step.
The right hand column of Table 4 is a pre-inversion conditioning sequence that includes steps for residual noise suppression, residual event timing (velocity) corrections, wavelet spectral equalization, residual amplitude balancing and data angle/offset muting of portions of the CRP gathers. A pre-inversion conditioning sequence is important prior to amplitude inversion because generic processing sequences (left column of Table 4) often do not adequately compensate for various acquisition and earth transmission effects.
Moreover for processing steps of both columns of Table 4, a user (typically a processing geophysicist) may be required to select parameters and parameter values for each of the processing steps to implement a particular noise or signal distortion correction. The selected values may critically alter the output amplitude characteristics of the data from that step. For instance, in a deconvolution step, the choice of operator type (gapped or whitening), the degree of time variation, the size of the autocorrelation design gate, the degree of spatial averaging and the length of the deconvolution operator can significantly alter the characteristics of the deconvolved data. The combination of selected parameters and parameter values from all the processing steps will determine the quality of the final amplitude inverted reflectivities.
Ultimately, the output data from the amplitude processing sequence should be optimized for amplitude inversion. Every step and every parameter in the processing sequence could be optimized by doing a global search for steps and parameter values on the whole data set. But for a 20 step sequence with 3 parameters per step and 10 possible values per parameter, the data set would have to be processed 1060 times to search the parameter value space—a totally impractical and cost prohibitive proposition. And beyond that, even if it were possible to generate 1060 output data sets, the problem of establishing criteria for selecting the optimized data set would remain. Consequently, because of the number of steps in a sequence, the order of the steps in the sequence, the number of critical parameters in each step and the amount of computing resources required to execute individual steps, the conventional practice is to test and select parameter values on a subset of the data one step at a time. An experienced interpreter or processor then judges which parameter values generated the highest quality output data for that step or set of steps by examining various quality control (QC) displays. The degree of spatial coherence in a stack section (sum of traces in CRP gathers), the coherence in a CMP/CRP gathers or the similarity of processed output to a theoretical response derived from well control are often the chosen QC displays for parameter selection. Because amplitude inversions are very sensitive to small lateral variations in amplitude between traces in CRP gathers, because outcomes from one step impact the effect of parameter values in later steps, and because typical QC products do not use outputs from an amplitude inversion, the parameter selection judgments for a single processing step can be suboptimal relative to the sequence goal of optimizing data for al amplitude inversion.
FIG. 9 shows single CMP and CRP gathers after a variety of processing steps have been performed including those steps included in Table 4. Note the significant changes in the S/N ratio and character of the signal (flattened events) as the data went sequentially through the indicated processing steps. FIGS. 10A & 10B compare an initial NMO corrected shot profile (FIG. 10A) to a final pre-inversion CRP gather showing that the S/N ratio (FIG. 10B), amplitude distribution and phase characteristics of events have been significantly altered. FIGS. 10A & 10B dramatically indicate how significant an influence a processing sequence has on the characteristics of the primary reflected signal that will be input to an amplitude inversion.
As those experienced in the art are aware, the seismic data on the left side (FIG. 10A) would generate noisy amplitude attributes. But even the data displayed on the right side (FIG. 10B) could generate noisy inverted amplitude attributes. This is because the data processing may not have adequately compensated for various residual acquisition and earth transmission effects. In fact when left uncorrected, small variations in an event's relative timing (velocity), amplitude and phase induced by such residual effects can lead to errors of 200% to 400% in an inverted attribute's (eg. gradient) RMS level. Even when such residual effects are quite large, they may be difficult to detect with typical processing QC methods that rely on a visual assessment of event continuity in CDP or CRP domains.
Following data processing, users must select various parameters for the amplitude inversion step itself (e.g. damping parameters and data weights). The choice of inversion constraint parameters often strongly impacts those reflectivities that are derived from differences in amplitude with offset. Users typically use two methods to quality control output from an amplitude inversion: examine common depth point (CDP) spatial coherence of inverted reflectivities or analyze how closely the inverted reflectivities compare to reflectivities calculated from existing well control. Neither of these methods indicates whether the selected inversion parameters and constraints are optimal nor do these criteria ensure that the choice of parameters in a zone distant from well control is optimal. In practice, inverted attributes are often subjected to further processing, statistical analysis and “calibration” to increase their signal to noise ratio. Such post inversion analysis often cannot compensate for inadequate processing without the introduction of additional uncertainty and error.
In summary and as described above, there are several shortcomings in conventional seismic processing and amplitude inversion. A first shortcoming is that an amplitude processing and inversion sequence may be less than optimal because of poor parameter selection, resulting in suboptimal signal to noise ratio of inverted reflectivities at and away from well control. A second shortcoming is the lack of criteria for QC (quality control) in the selection of processing and inversion parameters for zones of interest that are applicable at well control as well as distant from well control.
A third shortcoming in current amplitude inversion methods is the use of various constraints to supplement the information contained in processed seismic data. Preferably, the relative amplitude behavior of optimally processed CRP gathers is proportional to the angle dependent interface reflection coefficient and therefore requires minimal hard or soft constraint information in order to yield high signal to noise ratio (S/N) inverted attributes. When processed seismic data is of poor quality, the inversion algorithm itself may have to be strongly constrained (via damping and weighting parameters) to produce realistic values of inverted attributes. As constraint weights increase, the inverted attributes become more dependent on the constraints and less dependent on the processed seismic data. This brings into question the accuracy of inversion results away from locations where the constraints may not apply.
The present invention addresses these shortcomings.